Mon. Not. R. Astron. Soc. 000,[T]|5]() Printed 5 February 2008 (MN MgX style file v2.2) 



Algorithmic regularization with velocity-dependent forces 



Seppo Mikkola 1 * and David Merritt 2 f 

1 Tuorla Observatory, University of Turku, Vaisalantie 20, Piikkio, Finland 
2 Department of Physics, Rochester Institute of Technology, Rochester, NY 14623, USA 



5 February 2008 



ABSTRACT 

Algorithmic regularization uses a transformation of the equations of motion such that 
the leapfrog algorithm produces exact trajectories for two-body motion as well as reg- 
ular results in numerical integration of the motion of strongly interacting few-body 
systems. That algorithm alone is not sufficiently accurate and one must use the extrap- 
olation method for improved precision. This requires that the basic leapfrog algorithm 
be time-symmetric, which is not directly possible in the case of velocity-dependent 
forces, but is usually obtained with the help of the implicit midpoint method. Here we 
suggest an alternative explicit algorithmic regularization algorithm which can handle 
velocity-dependent forces. This is done with the help of a generalized midpoint method 
to obtain the required time symmetry, thus eliminating the need for the implicit mid- 
point method and allowing the use of extrapolation. 
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1 INTRODUCTION 

' In some Af-body problems one has velocity-dependent per- 
turbations. Examples are the relat ivistic terms, which are 
important in black hole dynamics llAarsetbJl200f) . or dissi- 
pative terms due to tidal friction or atmospheric friction 
in satellite orbits. The KS-reg u larization (e.g. basic KS: 

iKustaanheimo and Stiefel Il965t IStiefel and Scheifeld Il97ll 
and the CHAIN- method of iMi kkola a nd Aarsethlll993ft can 
easily handle any additional forces, however in multi-body 
regularization with the KS-transformation, large mass ra- 
tios cause problems. Therefore other regularization meth- 
ods -algorithmic regulari zations- such as the logarit h- 
mic Hamiltonian method (Mikko la and Tanikawalll999al lbl 

iPreto and Tremaindl 1999ft or the time-transformed leapfrog 
jMikkola and Aarsethl 12002ft must be considered. On the 
other hand, these methods, when combined with the ex- 
trapolation method iGragdll964 Il965l: iBulirsch and Stoerl 
1966) cannot easily include velocity-dependent forces, ex- 
cept with the help of the implicit midpoint method. Since 
implicit methods may be inefficient, there is motivation to 
study ways to make the integrations explicit, while at the 
same time utilizing the good properties of algorithmic reg- 
ularization. 

Algorithmic regularization is simpler than KS regular- 
ization and, what is most important, versions of it work for 



arbitrary mass ratios. This is especially important in sim- 
ulatio ns of black hole dynamics in galactic nuclei jMerrittl 

l20oeft . 

In this paper, we first introduce the problem using a 
perturbed two-body system as an example. Then we suggest 
a generalized midpoint method to be used as a tool to time- 
symmetrize any basic algorithm. Finally the generalization 
to the A'-body problem is briefly outlined. 



2 GENERALIZED ALGORITHMIC 
REGULARIZATION 

Here we discuss the formulation of the basic algorithms, the 
time-transformed leapfrogs, that are regular in two-body col- 
lisions. Then a generalized midpoint method, that can also 
be used with the Bulirsch-Stoer (BS) extrap olation method 
l|Gragdll964l . ll965l:lBulirsch and Stoerlll966ft . is introduced. 

2.1 The perturbed two-body problem 

We first consider the perturbed two-body problem with 
velocity-dependent forces. Let r and v be the position and 
velocity vectors respectively and m the mass of the two-body 
system and t the time. We may then write the equation of 
motion as 
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~ m ^3 + /( r > t >' U )> 
V. 
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(2) 
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This case is simple enough for a detailed discussion; general- 
ization to t he full iV-body problem will be straightforward. 

As iMikkola and Tanikawal Jl999albl) and 
iPreto and Tremaind (ll99St> demonstrated, there is a 
way to make the leapfrog algorithm exact for two-body 
orbits, and regular for two-body collisions in more compli- 
cated problems, if one introduces a time transformation. 
Here we concisely re-derive the algorithm and augment 
it to the case of a general (not necessarily Hamiltonian) 
perturbation. 

Let 



, m 1 2 

o = v 

r 2 



(3) 



be the binding (Kepler) energy of the two-body system. We 
have the energy equations 



1 2 . , rn 
-v + b — — 

2 r 



(4) 

b=-vf. (5) 

This allows the introduction of the two time transformations 

— 1 (c\ 

ds + () 



(7) 



dt r 
ds m ' 

which are equivalent along the solution trajectory. Using the 
first alternative Q to transform the equation of motion for 
the coordinates (t,r), one gets 



t = 



±w 2 + b' 
v 



and the second equation gives for b and v 

b' = -v g, 

i r , 

v = — j+ff, 



(8) 
(9) 

(10) 
(11) 



where primes indicate differentiation with respect to the new 
independent variable s and 



g=-f(r,t,v). 



(12) 



If the perturbation / (hence g) is independent of the velocity 
v, then the above equations allow the use of the leapfrog 
algorithm: 
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(13) 
(14) 

(15) 

(16) 
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the end of the step, and vi = (vq + v\)/2. If the per- 
turbation g — 0, then the motion is pure Kepler mo- 
tion and the leapfrog algo rithm produces an exact traje c- 
tory with only a time error l|Mikkola and Tanikawall999al lbl 
IPreto and TremainellT999l) . 

In the above equations, the symbol gi indicates 
g(ri). However, if g actually depends on the velocity 
too, then the leapfrog cannot be immediately formed. 
This problem (or rat h er an analogous one) was solved by 
Mikk ola and Aarsethl l)2002F) . using the implicit midpoint 
method, i.e. it was necessary to solve the equation 



Vi 



h—§- + hg I r i ,ti 



v + Vl 



(19) 



for vi. Often this solution is possible only by iteration which 
can be rather expensive if the perturbation is strong and 
complicated. This fact motivates a search for ways to find 
an alternative that is explicit, yet capable of utilizing the 
algorithmic regularization. This goal can be achieved with 
the help of the algorithm we next discuss. 

2.2 Generalized midpoint method 

Here we introduce a generalization to the well-known mod- 
ified midpoint method. In this algorithm, the basic approx- 
imation to advance the solution is not just the evaluation 
of the derivative at the midpoints, but any method to ap- 
proximate the solution. Thus the algorithmic regularization 
by the leapfrog can be used even when the additional force 
depends on velocities. That provides a regular basic algo- 
rithm, which is made suitable for the extrapolation method 
by means of the generalized midpoint method, as follows. 
Consider the differential equation 



z = f{z), z(0) = z . 
Splitting the above as 

* = f(y), 

V = f(x) 

with the initial values 

x = Vo =2(0), 

gives the leapfrog-like algorithm 

xi = Xo + ^f(y ), 
Vi = Vo + hf{xi), 



(20) 

(21) 
(22) 



Xi 



xi + -f(y 1 )- 



(23) 
(24) 

(25) 



However, this is nothing but another way to write the well- 
known modified midpoint method. 

A new interpretation of the above can be obtained by 
first rewriting it in the form 

xi = x + (+^/(Vo)) , (26) 



yi = Vo-[—2^ x \ 



Vi 



y± + (+^/( a; i)) . 



where the subscripts and 1 refer to the beginning and 



xi = xi - [--^f(yi) 



(27) 
(28) 
(29) 
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In 12611 the bracketed term is an (Euler-method) approx- 
imation to the increment of x over the time interval h/2 
with the initial value y Q , while in l|270 the initial value is 
xi ~ x(h/2) and the time interval is —h/2 Finally, this 
increment is added -with a minus sign- to y to obtain an 
approximation for y(h/2). In the remaining formulae I28H . 
I|29[l , the idea is the same but the roles of x and y have been 
changed. 

A generalization of this is now obvious. Let 



z(At) « z + d(z , At) 



(30) 



be an approximation to the solution of Eq. (121 H over a time 
interval At. In Euler's method, 

d(z ,At) = Atf(z ), (31) 

which gives the algorithm described in Eqs. 12611 - I29H . but 
in general, d could be obtained from any reasonable method 
for solving the differential equation 1201 . We thus choose a 
method and define 



d(z ,At) = z(At)-z , 



(32) 



where z(At) is the approximation for z(At) obtained with 
the chosen method. This generalized midpoint algorithm 
may be especially useful if one uses a special method that is 
well-suited to the particular problem at hand. 

One step in the generalized midpoint method can now 
be written 

d(Vo.+£), ( 33 ) 



x l 



Xq - 



y -d(xi 

Vi = yi+d(xi,+- t 

<*(»i, 



X i 



A(x,y,h) : x 



V 



x + d{y,+- 
y - d{x, -'- 



(34) 
(35) 
(36) 

(37) 
(38) 



2" 
2 j ' 

2 j ' 

"-) 
2 j ' 

or, if we define the mapping (or "subroutine" ) 

■-) 
2' 

■-) 
2 j ' 

we can write the algorithm with many (TV) steps as 

1. Set y — x; 

2. Repeat A(a;,y,/i) A(y,x,h) N times; (39) 

3. Accept x as the final result. 

Thus one simply calls the subroutine A alternately with 
arguments (x,y) and (y,x) such that the sequence is time- 
symmetric (starts and stops with x in Eg. 1391 . 

This basic algorithm has the correct symmetry - be- 
cause it was derived from a leapfrog-like treatment - such 
that the error in integration over a fixed time interval with 
different timesteps h can be written 



error = Aih 2 + A4/1 4 + . 



(40) 



and thus the Gragg-Bulirsch-Stoer extrapolation method 
can be used to obtain high accuracy. 

The great advantage of this generalized midpoint 
method is that the leapfrog with the implicit midpoint 
method can be replaced by a method that is not exactly 
time-symmetric. The computation of the quantity g 1 , when 




Figure 1. Number of perturbation evaluations per unit of time 
(dN/dt) in a two- black- hole system, with masses and speed of 
light mi = 0.9, m,2 =0.1, c = 20, integrated from the initial 
values ao = 1, eo = until the final merger of the two black holes. 
The i-axis is the semi-major axis in units of the Schwarzschild 
radius. The green curve is for the new method while red and 
blue illustrate two varieties of the implicit midpoint method (as 
described in the text). 



it depends on velocity, can be done in a straightforward way, 
e.g. by 



9i =9(ri,ti,vi] 



(41) 



where one may approximate vi either by v 1 m vo or prefer- 
ably by 



V L R!t)o- 



2 r\ 



(42) 



after which 



Vi = vo — ft-f- + hg(ri ,t 1 ,tii ) (43) 
7* ^ 222 

can be used instead of 1151 for ll9H . Here it is necessary to 
stress that only the increments of the variables from the 
algorithm I|13^ - ljl8|l are to be used as the quantities d in the 
algorithm jgfl-JScl. 



3 SOME EXPERIMENTS 

Using a simple perturbed two-body code, written accord- 
ing to the above theory, we carried out some experiments 
to compare the new alternative with the implicit midpoint 
method. 

Tests with an (initially) circular orbit of unit radius and 
with the perturbing (frictional) force / = — ev suggest that 
for very small e ( ^ 10 -6 ) the implicit midpoint method 
is faster, but for stronger perturbations, the new method is 
favorable. 

T ests with the relativistic PPN2.5 terms from ISoffell 
( 1989) are illustrated in Fig. Q Here the system was a two- 
body system with masses mi = 0.9, m,2 = 0.1, initial semi- 
major-axis ao = 1, initial eccentricity eo = and the ve- 
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Figure 2. Error evolution in the experiments. The measure of 
error, plotted as a function of time, is |(^« 2 + &)— — 1|. Colour 
coding is the same as in the previous figure. 



locity of light was set to c — 20. Due to the gravitational 
radiation term, the semi-major axis shrinks and the com- 
putational effort (dN/dt— number of perturbation evalua- 
tions per unit of time) increases. The figure illustrates the 
evolution of dN/dt (averages over 100 steps with BS extrap- 
olation) during the computation (until final merger of the 
bodies) for three different methods. The results are plotted 
as a function of the shrinking semi-major axis (measured in 
terms of the Schwarzschild radius for the combined mass). 
In these integrations the one-step relative error tolerance 
was set to 10" 13 and the errors, measured via the quan- 
tity -^(\v 2 + b) — 1, were ~ 10" 11 for the new method and 
the midpoint method with iteration to convergence (corre- 
sponding to the "green" and "red" experiments in Fig. 0. 
For the restricted iteration method the error was, however, 
about 10~ 9 suggesting that this method is not to be rec- 
ommended. The errors grew secularly, as can be seen from 
Fig. |3 the numbers given above refer to the values just at 
merger, i.e. when the two particles approach more closely 
than the sum of Schwarzschild radii. It may be seen that in 
all cases, the new method is somewhat more efficient. 

4 iV-BODY FORMULATION 

The generalization of the algorithm to the A-body prob- 
lem is simple in princ iple. One may use the l e apfrog algo- 
rithms introduced by i Mikkola and Tanikawal il999al lbh or 
Mikk ola and Aarsethl J2002I) and simply add the necessary 
velocity-dependent forces. A new formulation that effec- 
tively unifies the above cited works may be constructed 
as follows. Let T = (1/2) J^. mk,v\ be the kinetic energy, 
U — '^2 li< .m i vn,j\r i — Tj\~ x be the potential energy, and 
Q an (in principle) arbitrary function of of the coordinates, 
often 

n = 5^|r i -r i r 1 . (44) 

i<j 



Then one may define, in analogy with (|SJ and 0, the two 
time transformations 

t' = l/(aT + B) = l/(a[/ + /3fi + 7), (45) 

where a, j3 and 7 are adjustable constants. Since T = U + E, 
we have B — —aE + /3f2 + 7, which expression is used only 
for the initial value of B and later this quantity must be 
obtained by solving the differential equation 

B = -aJ2v k -f k + ^^-v k . (46) 

k k 

In the above, Vk, fk are the velocity and position of the body 
with mass rrik, correspondingly, and the forces additional to 
dU/dr k are denoted by f k . 

The equations of motion that can be used to construct 
the leapfrog that provides algorithmic regularization are, for 



time and coordinates respectively, 

t' = l/(aT + B), (47) 

r' k = t'vk (48) 
and for velocities and B 

r = l/(af/ + ^n+ 7 ), (49) 

v'k = r'(^- + f k )/m k , (50) 

B> = r>J2(-«fk + ^)-Vk. (51) 



Here the (possible) velocity dependence of the additional 
forces f k can be handled as in our two-body example 
above. However, to account for the (explicitely written) 
^-dep endence of B' one must follow IMikkola and Aarsethl 
(2002), i.e. first the v k are advanced and then the aver- 
age (wfe(O) + Vh(h))/2 is used to evaluate B' . Thus the 
leapfrog can be constructed in obvious analogy with the per- 
turbed two-body case. However, in A-body integrations, the 
roundoff can be a serious source of error and relative coor- 
dinates of close bodies must be used to reduce that effec t 
( Mik kola and TanikawUll999al : IMikkola and Aarsethll2002l) . 

Some additional remarks follow. 

(i) If one takes (u,f3, 7) = (1,0,0) then the method 
obtained is the logarithm ic Hamiltonian method 
l|Mikkola and Tanikawglll999ar i. 

(ii) If (a, j3, 7) = (0, 1,0) then we have the tim e trans- 
formed leapfrog (TTL) iMikkola and Aarsethir2002h . 

(iii) If (a, f3, 7) = (0, 0, 1) then the method is just the nor- 
mal basic leapfrog. 

(iv) If there are no velocity-dependent perturbations, then 
the normal leapfrog can be used and it is in fact faster. 
This is because our alternative algorithm then does some 
(unnecessary) calculations back and forth. 

(v) The question of which combination of the numbers 
(a, ft, 7) is best cannot be answered in general, but experi- 
mentation is necessary. For Af-body systems with very large 
mass ratios, however, it seems that one must have /3 / 0. 
which means a form of the TTL method. 

(vi) The experiments discussed in section |2] correspond to 
the alternative (i), i.e. (a, j3, 7) = (1, 0, 0). Note that for the 
case of only two bodies, there should be not much difference 
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between alternatives (i) an d (ii) since in this case they are 



between alternatives (1) ana (n) since m tms case tncy 
mathematically equivalent jMikkola and Aarsethll2002f) . 



5 CONCLUSION 

We have demonstrated that the generalized midpoint algo- 
rithm can be used to time-symmetrize the algorithmic regu- 
larization leapfrog even when the forces depend on velocities. 
This permits efficient use of the extrapolation method. For 
very small perturbations, the implicit midpoint method may 
still be better, and the new method can be recommended 
only when the velocity dependence of the forces is signifi- 
cant. Finally we note that the generalized midpoint method 
can be used with any special low order approximation to the 
differential equations under consideration. Thus it is not re- 
stricted to TV-body problems. 
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